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ABSTRACT 

It is generally accepted that the heating of gas in clusters of galaxies by active galac- 
tic nuclei (AGN) is a form of feedback. Feedback is required to ensure a long term, 
sustainable balance between heating and cooling. This work investigates the impact of 
proportional stochastic feedback on the energy balance in the intracluster medium. Us- 
ing a generalised analytical model for a cluster atmosphere, it is shown that an energy 
equilibrium can be reached exponentially quickly. Applying the tools of stochastic cal- 
culus it is demonstrated that the result is robust with regard to the model parameters, 
even though they affect the amount of variability in the system. 
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1 INTRODUCTION 

The cooling time of gas in the cores of many galaxy clusters 
is much shorter than the Hubble time. In the absence of 
heat sources this gas will cool and flow towards the centre of 
the cluster. However, high resolution X-ray spectroscopy has 
shown that the rate at which gas cools to low temperatures 
is much lower than i nitially expected (IPeterson et all l200ll ; 
iTamura et al.l 120011 ; IPeterson et al, I l2003h . suggesting that 
the gas is being reheated. 

Feedback from the active galactic nucleus (AGN) in 
the central galaxy of the cluster has been considered the 
most promising mechanism fo r the reheating of the cool- 
ing flow (see, e.g., review by iMcNamara fc Nulsenl 120071 , 
and references therein). Modelling of the heating process 
is difficult because the physical properties of the AGN- 
cluster interaction and AGN accretion are far from clear. 
It has been suggested that the cluster gas is heated by 
outflows, bubbles, sound waves, thermal conduction, turbu- 
lence, and (or) a combination of the above mentioned pro- 
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cesses (see, e.g.JHeinz et al.11 2006: Briiggen & Kaiser 2002; 



Ruszkowski ct al. 2004: Voigt et al. 2002; Fuiita et al. 2004; 



Scannapieco fc Briiggenll200ct ). Even though it is not clear 



how effective these individual processes are in raising tem- 
perature of the centre of the cluster, it is clear that the 
energy deposited by the AGN must be thermalised so that 
on average the heating and the cooling rates remain equal 
(or at least comparable). 

Currently, the theoretical description of AGN activity 
in cooling flow clusters is based primarily on results from 
numerical models. A typical numerical experiment would 
assume a (magneto-)hydrodynamical model for the intra- 
cluster plasma, and find a detailed numerical solution of the 
appropriate equations. In this framework the relative impor- 



tance of various specific physical processes can be evaluated 
and their role in the heating of cooling flows can be conse- 
quently assessed. It is prohibitively expensive, however, to 
cover a large parameter space in the simulations, and accu- 
mulation of statistics is often very difficult. While the sim- 
ulations shed light on some relevant physical phenomena, 
they are not well suited for answering questions regarding 
long-term stability, variability, and statistics of the popula- 
tion of clusters. 

In this article, an analytical model for the atmosphere 
of a cluster is constructed to study the stability and evo- 
lution of the energy balance driven by AGN feedback. The 
model behaviour is expressed in terms of stochastic differ- 
ential equations, which are widely used in studies of fluctu- 
ating phenomena, but can appear somewhat obscure at first 
glance. Consequently, the model is deliberately chosen to be 
very simple, so that the interpretation of the solutions of the 
model equations remains robust. This analytical approach 
is complementary to the more common detailed numerical 
simulations. It helps to probe otherwise difficult to reach pa- 
rameter space, and addresses some of the above mentioned 
questions regarding AGN feedback in clusters. Our approach 
to these questions in t he cur r ent w ork is closely related to 
the ideas developed by iPopel (|2007l ). The proposed answers 
are statistical in nature, which can be useful in selecting pa- 
rameters for simulations as well as interpreting their results. 

Methods of stochastic calculus have been used exten- 
sively in application to problems of statistical physics, chem- 
istry, computational biology and finance. In astrophysics 
they are mainly used in astroparticle research (see, e.g., 
lLitvinenkdl2009r ) and are quite suitable for investigating pa- 
rameter space in the context of phenomenological models, 
as demonstrated in the present work. 

The outline of the article is as follows. Section[2]presents 
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the main model assumptions and derives the equations for 
deterministic energy balance in the ICM. Section [3] intro- 
duces the concept of the stochastic mass deposition rate, 
discusses briefly the rules of the stochastic calculus, and de- 
rives stochastic differential equation for the heating of the 
ICM. Section [4] using a further simplification, solves the 
equation, derives formulae for the mean, the probability den- 
sity function of the heating process and its asymptotic form. 
Section [5] discusses possible observational implications and 
limitations of the presented model. 
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2 TOY MODEL 



Consider a relaxed cluster of galaxies. For simplicity the in- 
tracluster medium (ICM) in this model is assumed to be an 
unmagnetised, high temperature gas, confined by the grav- 
itational potential of the dark matter halo, which preserves 
its spherically symmetric distribution. 

The surface 53 of the sphere with radius r = r$ (see 
Fig. [T]) is considered to be the outer boundary of the ICM, 
beyond which the gas is no longer gravitationally bound to 
the cluster. The radial velocities and the densities in the re- 
gion r > r$ are assumed to be negligible, V3 = v(r > rs) ~ 0, 
9 — pij' > ?"3) « 0. The density of the ICM rises towards the 
centre of the cluster, and the high temperature ICM cools 
by emitting X-rays. The cooling of that gas results in the 
classical cooling flow - a spherically symmetric flow towards 
the centre of the cluster. The surface S2 of the sphere with 
radius r — T2 is located sufficiently far away from the cen- 
tre of the cluster so that the density and the temperature 
of the gas in the volume Voutcr remains unaffected by the 
AGN activity. Therefore, the absolute value of the velocity 
of the cooling flow on this surface «2 = v{r-i) is determined 
only by the cooling rate in the volume Voutcr, and remains 
approximately constant. The volume Virmer is heated with 
the integral rate H(t) 7^ 0, while it also cools with the in- 
tegral rate Li nncl (t) 7^ 0. In general, values of the velocity 
vi = v{ri,t), the density pi = p(ri,t), and the pressure 
Pi = p(ri,t) at Si are not constant because of the AGN 
activity. The region inside £1 is excluded from the model. 

It is straightforward to write the energy balance equa- 
tions for the volumes Vouter and Vinncr using the generic 
equation, 



dEy 

dt 



p— + pu + pip 



v.dS, 



(1) 



where Ey is the sum of kinetic, thermal and potential ener- 
gies of the gas inside the volume V, u> is the enthalpy of the 
gas, 4> is the gravitational potential, dS is the element of the 
surface with an external normal. If the velocity of the gas 
is subsonic almost everywhere, then v 2 /2 < uj, and we can 
neglect the contribution of the kinetic energy term in the 
equation fTJ. Energy balance for the volume Voutcr is then 
given by, 



(2) 



where M2 = 4nr 2 p2V2 is the constant mass deposition rate 
through S2. For the volume dinner, the energy balance equa- 



Figure 1. Division of the cluster into the inner and the outer 
regions. The following notations are used: H is the integral heat- 
ing rate, L is the integral cooling rate, n (i = 1,2,3) are the 
radii of the spheres, Si are the surfaces of the spheres, n; are 
the normal vectors, Vi nnor is the volume between the surfaces Si 
and S2, Voutcr is the volume between the surfaces S2 and 53. 
The orientation of the normal vectors ni and 113 is external to 
the considered volumes (note that ni is external to Vi nncr , which 
does not extend all the way to the centre) , whereas n 2 is external 
for Vi nn er and internal for Voutcr- 



tion can be written as 

H(t) - L inncl {t) = -(wi(t) +Vi)Mi(t) + {lu2 + ip2)M 2 
= — (wi(t) + tpi)Mi(i) + L 

outer 1 

(3) 

where Mi = 4-Kripi(t)vi{t) is the time dependent mass de- 
position rate through Si . Differentiation of both sides of the 
equation © gives, 



dt dt 



(u>i(t) + i>i) 



dMi 
dt 



duii / n 
-df M ^ 

L outcr )^^ (4) 
Mi dt 



= (#(*)- LfamerCt) 

where, in the last equation, the term —(u>i + ipi) was ex- 
pressed from the equation ((3|. It is reasonable to presume 
that H 3> Linncr, as the main contribution to the heating 
rate integral H is likely to come from a small volume of gas 
close to Si, in the vicinity of the AGN. Whereas the inte- 
gral cooling rate function has the entire volume of Vi nncr as 
a support, and, therefore, is likely to be less affected by the 
changes close to Si. Making this assumption results in the 
following differential equation for the heating of the ICM, 



dH 
~dt 



,„, , r , ., 1 dM duj •,, . 
\H{t)-L{t)]^-—-—M{t), 



(5) 



M dt dt 

where L(t) = L uter + dinner {t) , said subscripts '1' were 
dropped to simplify the notation. 
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If the heating matches the cooling in the inner region 
H — Li nn er = const, then equation ([5} reduces to, 



Lou 



1 AM 
KP~aT 



dP 



(6) 



which is easily integrable, 

[1 1 1 

Lo(t) — Lo(to) + Louter — : • (7) 

lM(t) M{to)\ 

In this case the enthalpy is inversely proportional to the rate 
of mass inflow. It is also straightforward to get this result 
directly from the equation ©. 

3 STOCHASTIC FEEDBACK 

AGN feedback process is thought to recycle the rest-mass 
energy of the infalling material during the course of its ac- 
cretion onto the central supermassive black hole. Details of 
the physics of accretion, methods of the energy transfer and 
thermalisation in the ICM are far from clear. Nevertheless, 
as the end product of a feedback cycle, the enthalpy of the 
ICM must rise and M must be reduced. The magnitude of 
the change of the mass flux, AM, is proportional to the mag- 
nitude of the AGN "response" , which in turn is likely to be 
proportional to the "input signal" M, 



AM oc M. 



(8) 



Such linear relation between the response and the input 
signal is found, e.g., in transient X-ray binaries, where the 
power of the outburst was found to be proportional to th e 
amount of the material in the disk (jShahbaz et al.|[l998f ). 
Although the physical processes responsible for this rela- 
tionship in X-ray binaries are likely to be quite different 
from the physics of AGN in clusters, it does suggest that 
the linear relation between response and the input may be 
a reasonable assumption. 

The proportionality coefficient in the equation ([8]| is 
unlikely to be a universal constant. In principle, magnitudes 
of the responses can vary given the same input M. Across 
a population of clusters such variation could be due to the 
differences in masses and spins of the central black holes. In 
an individual cluster, the constant of proportionality could 
also vary because of changes in the state of the accretion 
disk. 

In order to account for a range of possible responses, 
the value of AM (and, therefore, M) can be considered to 
be a random variable parameterised by time, i.e., a stochas- 
tic process. The equations of the toy model in combina- 
tion with a suitable distribution for M can then be used 
to infer properties of the AGN heating. The most unbiased 
choice of the distribution, which still reflects the property of 
the proportionality ([8]) is a uniform distribution of AM / M . 
The uniform distribution corresponds to the "white-noise" 
stochastic process ((t)At: (£(*)) = 0, (C(*)C(0) = S(t-t'). 
Using notations accepte d in the theory o f stochastic differen- 
tial equations (see, e.g.. lGardiner| [2005, for an introduction 
into the theory and applications of stochastic calculus) this 
choice of the distribution implies the following stochastic 
differential equation (SDE), 



where b is a constant, and W(t) is a Wiener process. 

The Wiener process is defined as a continuous time 
random walk (Brownian motion) in the limit of infinites- 
imally small step size. According to the rules of stochas- 
tic (also called ltd) calculus the equality ((t)At = AW(t) 
is interpreted symbolically as defining the integral relation 
W(t) — W(to) = f t ((t')At' , because strictly speaking the 
derivative of W(t) does not exist. The simplest way to find 
an analytical solution to a SDE is to make a suitable vari- 
able substitution. If the substitution transforms the term 
f(-)AW(t), where /(■) is a function, into cAW(t), where c is 
a constant, the equation can be easily integrated. According 
to the rules of stochastic calculus the change can be formally 
done by expanding the new variable to the second order, 
and using the following identities: (AW(t)) 2 = At, (At) 2 = 0, 
AW(t)At — 0, the first of which simply states that the vari- 
ance of Brownian motion is equal to the elapsed time. 

Expansion to second order of the function fi = InM 
relating to the stochastic process M is given by, 



(10) 



dw= idM ^(AM) 2 

M 2M 2 

= bAW(t) - ife 2 dt. 
The last equation can be straightforwardly integrated, 



H{t) = fi(t ) + b[W(t) - W(t )] - -b 2 (t - t ). (11) 

Changing the variable /j, back to M gives the analytical so- 
lution of the SDE ©, 

M(t) = M(to)e blW(t) - W{to)] - b2(t - to)/2 . (12) 

Using the fact that W(t) is a normally distributed random 
variable with zero mean and variance (W(t) 2 ) = t, it is easy 
to calculate the first two moments for M, 



(M(t)) = (M(t )), 
{M{t) 2 ) = (Af(to) 2 )e 62(t " to) . 



(13) 



AM 
~M 



bAW(t), 



(9) 



In other words, if AM /M is a white-noise process, then the 
mean of M remains constant, whereas the variance grows 
exponentially. 

Note that according to equation (|12|) the mass deposi- 
tion rate is non-negative, M > 0. The limiting value M = 
corresponds to the stable situation, when the heating rate 
equals the cooling rate exactly. According to the equation 
([21) if M = the heating rate is given by H = Li nncl + L ou t el . 
In the present model, however, the heating ability of AGN is 
explicitly limited to the inner part of the cluster. The max- 
imum heating rate is, therefore, H = Li nnel . If L ou t C r 7^ 
the mass deposition rate is always positive, M > 0, and the 
cold gas continuously accumulates in the cluster's centre. 

In order to understand how the stochastic behaviour of 
M influences the energy balance in the ICM, it is necessary 
to make a further assumption about the reaction of the ICM 
to the change in M. As was shown above, see Q , in the state 
of the stable deterministic equilibrium the enthalpy of the 
ICM, u, is inversely proportional to M. By presuming that 
this is also true in the case of the stochastically behaving 
M, the reaction of the enthalpy to the change in M can be 
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calculated using the rules of Ito calculus, 

1 1-2 
AiU = -Loutcr-r— dM + L ou t cr -t— (dM) 

M 2 M A 



L 



M 



(bdW - b 2 dt) . 



(14) 



Substitution of equations (|14|l and ([9]) into equation ([5]) 
yields the following SDE for the ICM heating, 

dH(t) = [H(t) - L inncv (t)]bdW(t) - b 2 Loutcrdt. (15) 



4 HEATING PROCESS 

Considering for simplicity that Li nnel ~ const, the heating 
SDE (|15|) can be rewritten with the dimensionless variables: 

h = H I ii nncr , A = ioutcr / -dinner i and T = 6 t, 



dh(r) = [h{r) - l]dW(r) - Adr. 



(16) 



To solve the equation (|16|l it is best to start by letting 
A = 0. This corresponds to the scenario when the AGN 
can heat the entire cluster. The change of the variable x = 
ln(l — h) leads to the following equation, 



da; = 



1 



.(dhf 



= dW(r) 



rdr, 



which can be easily integrated, giving the solution, 
h{r) = l-{l-ho)e W ^'\ 



(17) 



(18) 



where to = 0, ho — h(0), W(0) = 0. It is important to 
note, however, that although the change of variable guar- 
antees that h < 1 it also allows an unphysical situation 
when h < 0. To exclude solutions with the negative heat- 
ing it is necessary to treat h = as a boundary condition. 
Clusters that reach this boundary can be reflected at it (h 
turns positive), they can be absorbed at the boundary (h 
remains zero), or, possibly, a fraction of clusters can be re- 
flected and the remaining part gets absorbed. The reflective 
boundary seems to be a natural choice, since the case of the 
absorbing boundary requires permanent extinction of the 
active nucleus, while generally M 7^ 0, leading to continu- 
ous accumulation of the cold gas in the centre. According to 
the equation (|18[) the heating remains positive in the range 
W(t) < W = r/2 - ln(l - ho). Solving (TH} for W, and 
reflecting the function around the point W', gives the equa- 
tion for the reflected process, valid in the range W > W. 
The full solution of the heating SDE, with reflection at h — 
(W = W), is given by, 



h(r) 



1- (l-ho)e w{T) - T/2 , W(t)<W, 
1 _ -J-e- w W +T/! , W(t)>W. 



(19) 



Using the solution (JT9J) it is straightforward to calculate the 
mean value of the heating, 



(ft(r)> = 1 - ^erfc 



r/2 + ln(l - ho) 



2t 



h e T 



2(1 - ho) 



erfc 



3r/2 + ln(l - ho) 



2r 



(20) 



which shows that the equilibrium h — * 1 (H — > L) is ap- 
proached exponentially quickly as r — > 00, see Fig[2] Using 




Figure 2. Five sample heating curves (dotted lines) produced 
using equation (1201 ) with ho = 0.1, and the theoretical mean of 
the heating process (black line) as given by equation (|20fl . Note 
the change of the variability (dispersion) with time, compared to 
Fig. [3] 
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Figure 3. Two sample heating curves (dotted lines) produced 
using numerical solutions of equation II16I I with ho = 0.1, A = 0.2, 
and the mean of the heating process (black line) computed as an 
average of 10000 sample heating curves. Note that the amount 
of the variability (dispersion) remains constant in time, unlike in 
the case that A = 0, see Fig. [2] 



()19|) it is also possible to obtain analytical formulae for the 
variance and the autocorrelation function, but the expres- 
sions are excessively long and will not be reproduced here. 

The analysis of the solutions for the SDE can be com- 
plemented by considering the probability density function 
(PDF) of the processes. This is of interest for building a 
framework with which to interpret observations. For exam- 
ple, the most probable value of the observed heating rate 
is the modal value, where the PDF peaks. However, if the 
PDF is skewed, the mode might provide a poor estimate for 
the long term average heating rate, which is given by the 
mean of the distribution. 

The stochastic process described by the SDE (|16[1 has 
a PDF fhr(h, t), which satisfies the following Fokker-Plank 
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Figure 4. Evolution of the PDF for the heating process with 
A = 0, which started at r = as a delta function, S(h — 0.1). 
The curves correspond to the times r = 0.1, . . . , 5 with the step 
At = 0.1. 
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Figure 5. Asymptotic (r — > oo) shapes of the PDFs of the cluster 
heating according to the equations i !24[ i. If A 7^ the PDF has a 
maximum at 1 — A (in this case A = 0.2). For comparison both 
PDFs were normalised on the interval [0, 0.95]. 



equation (FPE), 

dfhr 
dr 



} df hT 1 d 2 
dh 2 dh 2 



[(h - lff hT ] 



(21) 



(see iGardinerl [2005, for a general derivation of the rela- 
tion between FPE and SDE). Substitution of the variable 
x — ln(l — h) and the corresponding change the PDF func- 
tion fh T = fx T e~ x (which preserves the probability measure) 
transforms the FPE (|21jl into the equation 



df XT 
dr 



0_ 

dx 



Ae" 



fxT 



+ 



1 d 2 f XT 

2 dx 2 : 



(22) 



which describes a diffusion process on the interval 1 x £ 
(—00, 0] in a medium with constant diffusion coefficient, 1/2, 
in presence of the inhomogeneous force field, 1/2 — \e~ x . The 
reflective boundary condition at x = can be satisfied auto- 
matically by extending the solution interval to the whole real 
axis, and reflecting the distribution at so that it becomes 
an even function. In the case where A = 0, the equation can 
be easily solved, 



Jxr — 



a x/2-r/8 



-(x-H) 2 /2r 



+ e 



-(z-K) 2 /27 



/«oe e/2 d£, 



(23) 

where /^o is the initial (r = 0) PDF. An example of how 
the PDF changes with time, starting from /50 oc <5(f — 0.1), 
is shown in Fig. [4] For large values of r, the PDF becomes 
a narrow peak at h = 1, in agreement with the SDE result, 
(h) -> 1. see (l20|l and Fig. [2] 

The connection between the FPE and the SDE also 
helps to explain the role of the (h— l)dW term in the origi- 
nal SDE If the term was simply dW the resulting FPE 
would be an equation of homogeneous diffusion. This would 
correspond to a case in which the feedback is not scaled with 
M. The limiting distribution in this case would be a uniform 
one. The h — 1 factor creates a force term in the equation, 
ensuring that, at all times, diffusion in the direction of h — 1 
is preferred over the opposite direction, no matter what kind 
of the initial distribution is assumed. 



1 Square bracket means that the end of the interval is closed, and 
the round bracket that that it is open. 



The FPE (|21|l can be rewritten in the form of a con- 
servation law, d T fhr = dhJhr, the function Jh T is called the 
probability flux function. In the stationary regime, the time 
derivative is zero, d T fh T = 0. It describes a statistically sta- 
ble scenario, which is presumably reached in the limit of 
large r. It follows that in this case Jht = const. Because of 
the reflective boundary at h — the fluxes on the left and 
right of this point must have opposite signs Jot = — Jot = 0, 
and therefore in the stationary case the probability flux is 
zero. Jh T = is an ordinary differential equation, which can 
be readily solved, 



fh oc 



-2A/(l-/i) 



A = 0, 

A / 0. 



(24) 



In the case A 7^ the normalised distribution is given by, 



fh 



2A 



(1 - h) 2 



(25) 



In the case A = the function fh is not normalisable on 
[0, 1]. However, the solutions (|24[> are still useful for under- 
standing the qualitative distinction between the two cases. 
If A / the PDF has a peak at h = 1 — A, see Fig. \E\ sug- 
gesting that clusters in this case would appear to be under- 
heated, on average. The change in the character of the PDF 
is due to the inhomogeneous force 1/2 — A/(l — h) in the 
diffusion equation (1221) . Along with the constant term 1/2, 
which forces diffusion in the direction of h = 1, there is a 
counteracting term, — A/(l — h), which reverses the diffusion 
on the interval 1 — 2A < h < 1. 

A physical interpretation of this property becomes 
transparent in the context of the toy model. The ability 
of the AGN to heat the cluster is limited to the volume 
Vinner (with X-ray luminosity Li nncr or 1). The additional 
mass deposition from the volume K> u ter (with the luminos- 
ity Louter or A) can not be counterbalanced by the AGN 
heating. This results in a shift of the balance from the exact 
match, h = 1 + A, to the lower end with the most likely value 
h = 1 — A (H — Linnet — Louter). The shape of the PDF is 
asymmetric and varies with A, see Fig. [6] Also the mode or 
the most probable heating rate is different from the mean 
or the average heating rate in this the distribution 

is skewed. 

It is possible to find an analytical solution to the SDE 
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Figure 6. Normalised PDFs in the case A ^ as given by 
the equation ( 12.">[) . The curves correspond to the values of A = 
0.02, . . . , 0.92 with the step 0.1. 

(|16|) in the case A / 0, 

h(r) = l-e(r)^ho + x£ , (26) 

where e(r) = e ( T )~ T / 2 is known as the exponential mar- 
tingale. This solution must be completed with the reflected 
process, so that h(r) > 0. However, because the analyti- 
cal quadrature of the integral term is unknown (see, e.g., 
iGoodman fc Kimll2006l . for details), this can not be done as 
easily as in the case A = 0. 

The long-term behaviour of the process in the case 
A / was already studied using the approach based on 
the FPE (|22p . The sample heating curves can be calculated 
using numerical solutions of the SDE Q5) . In Fig.[3]two sam- 
ple curves are plotted for processes with A = 0.2, starting 
from ha = 0.1. The mean of the heating rapidly approaches 
the value (h) « 0.55 (note that h — 0.8 is the most probable 
value) while contrary to the A = case the standard devi- 
ation remains finite, and the sample curves exhibit a lot of 
variation. 

Note that the balance of heating and cooling in the case 
when A 7^ is given by h — 1 + A (not h = 1) and therefore it 
is never reached. This is a direct consequence of the imposed 
limit on the AGN "fuel supply" (i.e., its heating capability, 
see section [2} to the region Vinner- The mass deposited from 
the Voutcr is not a part of the AGN feedback cycle, and 
accumulates in the centre. In practice this sets a limit to the 
value of A, which could be determined from future cluster 
surveys. 



5 DISCUSSION 

In the context of the model developed here, it is not sur- 
prising that most observed clusters are found to have a very 
tight correspondence between the heating p ower and the ra- 
diative cooling rate (|Dunn fc Fabian 20Q3)- While it would 
be naive to equate the heating rate in the current model 
with the power input measured from the size of the cavities 
in the clusters, it is reasonable to think that the two energy 
rates are positively correlated. 

In the presented model the self-tuning of the energy bal- 
ance is rapid, exponentially quick, in fact tx e T (see equation 



(|20[) ~). It works even in the case when the AGN is able to heat 
only a fraction of the cluster volume (see Fig. [3}. This shows 
that the correspondence of the cooling power and the AGN 
heating seen in the observations may indeed be a general 
rule, which is characteristic of the feedback process. Since 
the underlying PDF of the heating can be skewed (see Fig. [5] 
and equation (|25|) ~) it is worth noting that the heating rates 
inferred from the observations most frequently must lie near 
the mode of the distribution, and can not be used directly to 
find the mean or the average heating rate. The true picture 
can only be uncovered in a complete survey. 

As was shown above, the time dependent PDF becomes 
very close to the static solution for times r > 10, see Figs. 
[2] and [3] The physical time is proportional to r with the 
proportionality constant b~ 2 (see the normalisation of vari- 
ables just before the equation 1)160 ). The physical time limit 
is, therefore, t > 10fo~ 2 , which puts the lower limit on the 
model parameter b. If the physical time is of the order of the 
Hubble time the limit on b is given by b > 3 x 10" 5 yr _1/2 . 
In principle, this value can be verified observationally using 
the relation between the variance of M at different times, 
as given by the equation (|13|l . Given sufficient statistics for 
the mass deposition rates across the population of clusters 
in different redshift bins, the prediction of the exponential 
growth of the {M(t) 2 } can be tested. The values of the model 
parameters b and A can also be found. It is unlikely that, 
over the entire lifetime of a cluster, the real distribution of 
the heating is very close to the A = case. Since in this case 
the AGN would, on average, prevent any cooling flow from 
developing and also prevent any star formation, which, in 
reality, is ob served to be enhanc ed in some brightest clus- 
ter galaxies (|Raffertv et al.ll2006l ). This leaves the single pa- 
rameter distribution ()25p as a plausible heating PDF for the 
population of "mature" (of age t > 10) clusters. 

The current model sketches just one possible scenario 
for the feedback. In order to make the model robust and sim- 
ple, many assumptions were made: the insignificance of the 
magnetic field and the kinetic energy, rapid thermalisation 
of the AGN energy, and the inverse proportionality between 
the enthalpy and the mass flow. They have to be tested and 
verified in the future using numerical simulations and anal- 
ysis of the observations. In addition, it is quite possible, for 
example, that the re can be a delay in the AGN response to 
the cooling flow l|Popell2007h . During this delay an excess 
of cold material could accumulate around the AGN. This 
would be equivalent to effectively growing the cluster in size 
from the point of view of the feedback loop. AGN heating 
of such a larger cluster could then result in the overheating 
of the ICM, and possibly act as a destabilising factor. 

Despite the limitations of the model, this work demon- 
strates that the employment of stochastic calculus appears 
to be particularly well suited for the analysis of AGN feed- 
back. The combination of Ito calculus and the Fokker-Plank 
equations provides a simple, yet powerful, way of investi- 
gating behaviour of the time-dependent variability of AGN 
feedback in clusters. 
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